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Abstract 

In this paper, we construct a semi-implicit finite difference method for the time dependent 
Poisson-Nernst-Planck system. Although the Poisson-Nernst-Planck system is a nonlinear sys¬ 
tem, the numerical method presented in this paper only needs to solve a linear system at each 
time step, which can be done very efficiently. The rigorous proof for the mass conservation and 
electric potential energy decay are shown. Moreover, mesh refinement analysis shows that the 
method is second order convergent in space and first order convergent in time. Finally we point 
out that our method can be easily extended to the case of multi-ions. 
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1. Introduction 


The classical unsteady dimensionless drift-diffusion system which describes the evolution of posi¬ 
tive and negative charged particles p , n, and electric potential is given as follows [1] 

f p t = V • (Vp +pV</>), in H T := (0,T] x H, 

< nt = V • (Vn — nV</>), in Hr, (1) 

I — A</> = p — n, in Hr, 

where H is a bounded domain, [0, T] is the time interval, the length scale is chosen as the Debye 
length, and the time scale is chosen as the diffusive time scale. Note that the Debye length is much 
smaller than the physical characteristic length in most cases. Thus, if the length scale is chosen as 

the physical characteristic length scale in these cases, there will be a small parameter in front of 

the electric potential term in the Poisson equation (1), which will result in a singular perturbation 
problem [2-5]. However, in this paper, we only consider the case that the characteristic length scale 
is the same order of the Debye length, which also has some applications. For example, inside the 
ion channel in the cell membrane, the characteristic length scale of the ion channel is the same 
order of the Debye length. In this situation, we could choose the Debye length as the length scale 
so that the dimensionless system (1) is meaningful. The above system is called Poisson-Nernst- 
Planck system, which was first formulated by W. Nernst and M. Planck to describe the potential 
difference in a galvanic cell. The system has lots of applications in electrochemistry [6], biology [7] 
and semiconductors [8-10]. Based on the analytical derivations for the energy and entropy laws, 
Schmuck theoretically proved the existence, and in some cases uniqueness of the weak solution for 
the more general context of the Navier-Stokes-Nernst-Planck-Poisson system [11], In this paper, we 
focus on the numerical solution of the Poisson-Nernst-Planck system (1). 

When numerically solving the PDEs, to keep the original physical feature is greatly important 
in constructing numerical schemes for different physical problems. For example, one successful and 
active research is to construct structure-preserving scheme for the ODE systems (see [12] and refer¬ 
ences therein). Only schemes that are carefully designed can preserve mass and energy conservative 
properties. For example, finite difference methods were developed for solving the Euler Equations 
and Burgers equations to preserve the discrete energy dynamics [13]. In [14], authors used a finite 
volume method for the shallow water equations which conserves the mass, momentum and energy of 
the system. A finite difference method was presented in [15] for solving the nonlinear Klcin-Gordon 
equation which preserves the total energy. Qiao et al. [16] showed an unconditionally stable finite 
difference scheme for the dynamics of the molecular beam epitaxy, where the scheme preserves the 
energy decay rate exactly at discrete level. In [17], authors developed a general method of discretizing 
PDEs that preserving the energy using the average vector field method. Chiu et al. [18] developed a 
general mesh-free scheme for solving PDEs that can preserve the energy at discrete level. Chang et 
al. [19] discussed conservative and nonconservative properties of eight finite difference schemes for 
solving the generalized nonlinear Schrodinger equation. Chen et al. [20,21] proposed energy conser¬ 
vative finite difference schemes for solving the 2D and 3D Maxwell equations, respectively. 

In the past several decades, there appears a wide range of literature on numerical methods for 
the Poisson-Nernst-Planck system, including finite difference method, finite element method, and 
finite volume method, see [22-31] and references therein. Here we just mention some of the recent 
work, Bessemoulin-Chatard [29] gave a conservative finite volume method for solving the drift- 
diffusion equation, where the entropy inequality is preserved. Flavcll et al. [30] and Liu et al. [31] 
constructed two different conservative finite difference methods which satisfy the mass preserving, 
ion concentration positivity as well as total free energy dissipation numerically, where the total free 
energy is related to both electric potential and ion concentration, which is called entropy in [1], And 
Prohl et al. [1] presented two different finite element methods which satisfy electric potential energy 
decay and entropy decay properties, respectively. Now we briefly illustrate the first part of the work 
in [1] as follows: under the following initial conditions and zero Neumann boundary conditions 
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p( 0, x) = Po(x) > 0, n( 0, x) = no(x) > 0, in Q, 


(2) 


^ = ^ = ^=0, ondn T :=(0,T}xdn. (3) 

On on On 

It is well known [9] that the non-negative p, n is conserved in Qt, and the system (3) satisfies mass 
conservation, that is, for any t £ (0, T], 


M p = / po(x)dx = / p(t,x)dx, M n = / no(x)dx = / n(t,x)da 


( 4 ) 


where M p , M n are two positive constants which must be the same, since from (1) and (3) we have 


M p — M n 


/ (p — n)dx 
Jn 


f |ti S = 0 . 

Ian on 


And it is also shown in [9] and [11] that the system satisfies the following energy law 

E(t) + f f ((p — n ) 2 + (p + n) \Vcj)\ 2 )dxdt = E(0), (5) 

Jo Jn 

where E(t) = \ f n | V<p\ 2 dx is the electric potential energy. The above energy law (5) can be rewritten 
as 

dj? r 

= - y ((P - nf + {p + rOlWI 2 )dx. (6) 

Prohl et al. [1] proposed a finite element method which can preserve the mass conservation (4), 
ion concentration positivity and electric potential energy decay (5) in [1], However, the scheme in [1] 
is fully implicit, one has to solve a nonlinear system at each time step. In [1], a fixed point iteration 
method is used to solve the nonlinear system at each time step in order to get the rigorous physical 
quantities preserving results. In this paper, we present a simple semi-implicit finite difference method 
for the Poisson-Nernst-Planck system. For the new scheme, the unknown variables at next time step 
form a linear system which can be solved efficiently, no iteration is needed. Furthermore, the new 
scheme preserves mass conservation and electric potential energy identity numerically. Numerical 
results confirm the above properties. Mesh refinement analysis shows that the method is second 
order convergent in space and first order convergent in time. Finallly, we point out that our method 
can be extended to the case of multi-ions without any difficulty. 

The rest of the paper is organized as follows, section 2 gives the detailed numerical scheme and 
its properties, section 3 discusses the extension of the method for the case of multi-ions, section 4 
shows the numerical results, and conclusions and discussions are given in the final section. 


2. Numerical method 

In this section, we will develop a finite difference method which can guarantee the mass conser¬ 
vation (4) and energy decay (6) numerically. 

Although the method presented in the following can be extended into three dimension without 
any difficulty, we only give a detailed description when is a two dimensional rectangular domain, 
i.e. f l = [a, 6] x [c, d]. Let N x . N y be positive integers, the domain is uniformly partitioned with 
Ax = , Ay = and variables are stored at each cell center as follows 

^h = {(%j,Vk)\xj = a+ (j - i)A x,y k =c+(k- ^)Ay, 1 < j < N x , 1 < k < N y }. 

And time step is denoted by At. 

For a given two-dimensional grid function fj j., we define the following difference operators: 
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t -7 f t Jj+ Ijfc Jj,k fj,k+ 1 fj,k •, 

* hjj,k = (-—-, -—- ), 


Ahfj,k — 


As ’ Ay 
fj+l,k 2 fj t k ~b fj—l,k fj,k+l ~ 2 fjk + fj,k—l 


(As) 


X f _ fj+ 1/2, fc fj—l/2,k x f 

o x ;j, k - Ax , <Vi,fc- 


(Ay) 2 

fj,k+ 1/2 — fj,k- 1/2 


Ay 


d*hk = 


fj,k — r _ fj,k fj,k— 1 


As 


. d y fj, k = 


Ay 


where 


r _ fj,k + fj+l,k „ _ /j,fc + fj,k +1 

Jj-\-l/2,k 2 5 «/ 5 ^ ~I - 1/2 2 

The discrete L 2 inner product and the discrete L 2 norm are defined as 

< f,9>h= EE fj,k9j,kAxAy, || / ||h=< fj> h - 

3= 1 fc=l 


We also define 


iV* AT b 


< V/,Vy>^= EE Vfcjj.fc ■ S7 h gj,kAxAy, || V/ ||£=< V/, V/ > ?l . 

j=i fc=i 


2.1. Description of the method 


For zero Neumann boundary conditions (3), we define the values on center of the fictitious cells 
outside the boundary as follows 


PO,k = Pl,k,Pj,0 = Pj,l,PN x + l,k = PN x ,k,Pj,N y +l =Pj,N y , 

^0 ,/e n ltk ,n jt o rij^\,TlN x J r \^k ^Nx^ki^j^Ny+l ^j^Nyt ( 7 ) 

00 ,k = 01,A;, 0+0 = 0J,1,0AT X + 1,A = 0AT a ,,fe, 03,AT t/ +l = <fij,N y - 

where values on the the fictitious cells are denoted by the subscript with 0, N x + 1, and N y + 1. 
Our scheme for the Poisson-Nernst-Planck system (1) is as follows 


n m+1 

Pj,k 


Pj-.k 


At 


= A h p 


■t+l/2 


j,k 


+ 


(5 X (p™ k S X <f ™ + l/2 ) + Syip^Sy^ 112 
for j = 1 • • • N x , fc = 1 • • • ./Vy, to > 0 


( 8 ) 


where 


m+l _ m 

^ J ’ fc = A h n™+ 1/2 - (<5,K fc ^0™ fc +1/2 ) + WA€* +1/2 )) 


At 


for j = 1 ■ • • /V x , fc = 1 • • • iVy, to > 0, 


-A/i0”\ +1 = p”\ +1 - <+\ for j = 1 ■ • • N x , k = 1 • • • N y , m > -1, 


,m+1/2 

07,fc 


dm _i_ j.m+1 „ra _i_ m+l 

Pj,k _m+l/2 _ "T" ^’,/c _m+l/2 

‘Pj,k o 


^k + ^k 1 


(9) 

( 10 ) 


2 2 ■ 7 ’ K 2 
We should mention that the discrete equation for electric potential (10) is used for m = —1, since 
there is no initial conditions for 0, see (2). 
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2.2. Implementation of the finite difference method 


The method (8)-(10) is a semi-implicit method, which can be implemented efficiently. Assume 
that the quantities p , n, (j> are known at the previous time step m, and rewrite equations (8)-(10) in 
matrix and vector form as follows 



F IF <j) m + 1 4- 

_ )p m+l = (_ + _ )p m + A{pm) - +- 


( 11 ) 



—)N m+1 = 


'At 


-)A" 


A(N m ) 


$ m +1 -|- 
2 


( 12 ) 


F $ m+1 = N m+1 - P m+1 , (13) 

where I is the identity matrix, F is the matrix form of discrete Laplacian operator A^, A(P m ) is 
the coefficient matrix obtained from the second term of (8), which is a linear operator for P m , and 
P m+1 , N m+1 , $ m+1 are vector form of p, n, <f> at time step m + 1. 

As we can see that (11)-(13) is actually a linear system of unknowns P m+1 , iV m+1 and <f> m+1 , we 
can eliminate P m+1 and N m+1 in (13) using (11) and (12), this yields 


((^1 - F)F + A(P m + N m ))<f> m+1 = (^ J + F)(N m - P m ) - A(P m + N m )<f> m . (14) 

Once 4> m+1 is obtained, we can obtain P m+1 and N m+1 by solving (11) and (12), respectively. 

As mentioned in the introduction, p, n are non-negative in the entire time interval [0,T] as 
theoretically shown in [9]. In addition, through lots of numerical tests, it is found that the numerical 
solutions of p m : N m from the scheme (8)-(10) are also non-negative. If we make an assumption 
that the finite difference approximation of the derivative of <f> m is uniformly bounded in the entire 
computational time as done in [30] (see equation (37) in [30]), then we can prove that the numerical 
solutions of p m ; ]V m are also non-negative under some restriction of the time step size and mesh 
size followed by a similar proof of [30]. However, in general, it is not suitable to make such a prior 
assumption for the numerical solution which involves the quantities in the next time step to prove the 
properties of the numerical solutions. Thus, it will be a very hard task to prove the non-negativity 
of p rn i N m without the prior assumption for the finite difference approximation of the derivative 
of 4> m . But a lot of numerical tests show that the numerical scheme (8)-(10) produces non-negative 
P m , N m . Thus, in the following discussion, we simply assume that the numerical solutions P m , N m 
are non-negative. 

Theorem 2.1 Within the numerical solution of <f> m+1 up to a constant, the numerical solutions 
pm+i jjV m+i ;$ m+i 0 f fg)-(lO) are unique. 

Proof From (8)-(10), we can see that F, A(P m ) and A(N m ) are all symmetric banded and have 
the same matrix element structure. Moreover, for each of these three matrices, the diagonal elements 
are negative while the sum of each row is zero, we can get all eigenvalues of each matrix are less 
than or equal to zero through Gerschgorin Circle Theorem [32]. Thus all these three matrices are all 
negative semi-definite. Indeed, each of these three matrices has exactly one zero eigenvalue and all 
other negative eigenvalues. Furthermore, we can get ^ — F is positive definite. Since F is negative 
semi-definite, — F is positive definite, and F can be exchanged with ^ — F, we have (^ — F)F 
is negative semi-definite. And from above, we have A(P m + N m ) is negative semi-definite. Thus 
(■^ — F)F + A(P m + N m ), the coefficient matrix of (14), is negative semi-definite. And it is easy to 
see that (^ — ^)F + A(P m +N rn ) has exactly one zero eigenvalue and all other negative eigenvalues. 
Therefore within the numerical solution of $ m+1 up to a constant, the numerical solution of (10) is 
unique. 

Once we get $ m+1 , we can get P m+1 and N m+l by solving (11) and (12), respectively. Since 
-gj — is positive definite, the numerical solutions of P m+1 and N m+1 are unique. This completes 
the proof. □ 
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Since the coefficient matrix of (14) is negative semi-definite, symmetric and banded while the 
coefficient matrix of (11) and (12) is positive definite, symmetric and banded, these linear systems 
can be numerically solved very efficiently. 

In numerical computation using (8)-(10), we set cf> m+1 to be zero at one boundary point at each 
time step for (10) so that the (f> m+1 is uniquely determined. 


2.3. Main properties of the numerical scheme 


Theorem 2.2 For the solutions of (8)-(10), the discrete form of mass conservation (4) holds, that 
is, for any m > 0, 

N X Ny N X Ny 

E E^ lA * A y = ( 15 ) 

j= 1 k— 1 j= 1 k—1 

and 

N X Ny Nx Ny 

EE n r* +lAiA » = EE n™ k AxAy. (16) 

j= 1 k— 1 j=l k= 1 

Proof Multiplying At Ax Ay to both sides of (8) and summing for j = 1 • • • N x , k = 1 • • • N y , and 
applying the boundary conditions (7), we get 

N x Ny N x Ny 

E E -EE r?* A * A v 

j= 1 k= 1 j= 1 k=l 

N X Ny N x Ny 

= At Ax Ay ^2 AhpjY 1 ' 2 + At Ax Ay ^ ^ 4fa™A<Cfc +1/2 ) 

j= 1 fc=l J —1 fc=l 

+ At Ax Ay EEtKAC")' < 17 ) 

l=i fc=i 


From (7), the first term of the right hand side of (17) is zero obviously. The second term is 


N X Ny 

At Ax Ay £ J! S *(Pj:Mu L/Z ) 
j=i fc=i 

At Ay Ny Nx 


Ax 
At Ay 


+ 1 / 2 ' 
k 

m ( i m +h m /i m +5 /.m +3 

Pj+ %, k {<P j+ l*k 2 ) I ;fc W>pfc 

fc=l 1=1 

1V„ 

PjV. + i,fe(0iV,+l,fc ^ $N x ,k ) “Pi fcfal.fc - <A 0 ,fc 


EE( 

fc=i i= 

Ny 

Ax ^( : 


fc=l 


(18) 


_ ^ ^ TYZ | — TTt | — TYl~ j - “ TTt | ~ 

For zero Neuman boundary conditions, we have <f> 1 k 2 = 4> 0 k 2 and + 2 , = <p N k . Thus the 

second term is zero. Similarly, the third term is also zero. Therefore, the mass conservation identity 
for p is proved. 

Similar proof can be used for n. This completes the proof. □ 

Theorem 2.3 For the solutions of (8)-(10), the discrete form of energy identity (6) holds, that is, 
for any m > 0, 

pra+l _ pm 2 _ 2 _ 

-~" i/o vPL+<d*r +1/2 - v / pf+^id y v m+1/2 


At 


p m+ 1/2 _ n m+l/2 


(19) 
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where E m = 5 || V0 m || 2 , p™ and p R denotes the average value of p m in x and y direction, i.e. 

n m 1 m m , m 

t„m\ Pj-l,k'Pj,k t„m\ Ppk-l'Pj,k 
\PL )j,k - - 2 -’ (P/lW.fc - - 2 -’ 

and similar for nff and nff. 

Proof. Equation (10) of time level m minus equation (10) of time level m — 1 gives, 

-A ,^ 1 - <f% k ) = (p ™ +1 - p% k ) (n ™ +1 - n&). ( 20 ) 

Multiplying above equation with AxAy to both sides and summing for j = 1 • ■ • 7V X , k = 

1 ■ • • lV y , and applying the boundary conditions (7), we get 


< 


p m +1 _ pin 


At 


„m+l _ „m 

^ +1/2 >.-< A , ,^ +1/2 > h 


At 


N x N, 


h EE A ^(C fc +1 - C fc )(C fc +1 + C fc ) A ^ A ?/ 


2At 


1 

2At 


l=i fc =1 
IV* JVy 

E E V *(C * +1 - Cfc) ■ V *(C * +1 + C*) A * A y 

1=1 *=i 

IV* JVy 


= ^EE iv*c* + 1 i 2 - iv,c*i 2a ^ 


1 

2A t 

3=1 k =1 
£™-+l _ pjm 


At 

Substituting (8) and (9) into above equation, we get 
E m+ 1 — E m 


( 21 ) 


At x A h p m+1 / 2 + (4 (p m M m+1/2 ) + S y (p m <5 y <T +1/2 )) , <F +1/2 ) h 

- (A^n m+1 / 2 - (4 (n m 5 x cT +1/2 )+5y (n m 6 y <T +1/2 )) , <T +1/2 

= ((p m +v 2 , A h d> m+i ' 2 ) h - 1| vmv 


jm+ 1/2 


^n m+1 ' /2 ,A /i ^ m+1 ' /2 |) 


mahim+l /2 






\/^ y V 


frf™r) h rh m+1 / 2 

n R u cp 


= (p m+1/2 -n m+1/2 ,A^ m+1 / 2 \ - Vp™ + <9^ m+1 / 2 - VPR+riRdyt m+1/2 

( 22 ) 

Here we have used 

^A h p m+1 / 2 , <(> m+1/2 } ;i = - (V^p m+1 / 2 , V fc ^ ra+1 / 2 ) h = (p m+1/2 , A^ m+1 / 2 )^ , 

^A^n^ 1 / 2 , </> m+1/2 ) ?i = - (V h n m+1 / 2 , S/ h (j) m+1 / 2 ^ h = (n m+1 / 2 , A h <l> m+1/2 ) , 

and 


(23) 

(24) 


4 

(p m 5 x cT +1/2 ) 

,ci m+ i/2\ =- 

/ h 

VKd h *r +1 ' 2 

2 

7 

h 

(25) 


(p m 5 v cT +l/2 ) 

,m+l/2\ = _ 

/ /i 

v^<T +1/2 

2 

5 

h 

(26) 

s x 

(n m M m+1/2 ) 

^m+l/2\ = _ 

/ 7i 

^a x V m+1/2 

2 

? 

h 

(27) 
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(Sy ,^ m+l/2 } h = - 


l R u y 1 


i+l/2 


(28) 


Equations (23)-(28) can be easily checked when applying the boundary conditions (7). 

Equation (10) of time level m plus equation (10) of time level m — 1 gives, 

-A + Cfc) = G& +1 + PTk) K+ 1 + n™ k ). (29) 

Substituting (29) into (22), we get equation (19). This completes the proof of Theorem 2.3. 

3. Extending the method to the Poisson-Nernst-Planck system with multi-ions 


In this section, we shall extend the above method to the case of multi-ions. The model equations 
are as follows, 

(cj)t = V • (Vcj + ZidS7(j>), in O, (30) 

—A (f) = '^^ z iCi, in U, (31) 

i 

where d is a ion with valence Zi. Using a similar method as in the introduction part, it is easy 
to check the above system satisfy the following energy and mass identities under zero Neumann 
boundary conditions, 

-TT = - I (|£^r + C£,*tci)m 2 )dx, 

J £1 i A 


d /* 

— / d(t,x)dx = 0. 
dt Jq 

The scheme for the above Poisson-Nernst-Planck system is as follows 


+) 


At 


m 

— = Ma)Tk 1/2 + * (E (( Ci )r fe 4C fc +1/2 ) + Sy (te)MC* +1/2 )) > 




(32) 

(33) 


(34) 

(35) 


where 


,m+l/2 _ Vj,k 
P j,k ~ 




( \m+l /2 _ 

\ c i)j,k ~ 

The matrix and vector form of scheme (34)-(35), is as follows, 


(Ci) 


m+1 


(Ci) 


life 


IF , IF + <J>™ 

(At - 2 '^ At + 2 >^ + -7—’ 


F ^m+l = Zi C\ 


m+1 


(36) 

(37) 


(36) and (37) are linear system of C™' +1 and <f> m+1 , which can be solved efficiently, since all these 
matrices are symmetric and banded. 

Using similar techniques as in Theorem 2.2 and 2.3, we could also prove that the above numerical 
scheme satisfies the mass conservation and energy decay properties, which is stated in the following 
theorem. 


Theorem 3.1 For the solutions of (3f)-(35), the discrete form of mass conservation (33) and 
energy identity (32) holds, that is, for any m> 0, 

N x N y N x N y 

E X>)£ * lAxA y = EE^-)u lrA «- 

j=1 k= 1 j= 1 k=1 


(38) 










Moreover, the following energy identity is preserved: 


E m +1 — 


2 

yt^ c ^ m+1/2 

At 

2 

h 

2 

- 

^ z ? c ?R d & m+1/2 

5 

h 



(39) 


where E m = ||V<(> m ||)(, c™ L and denotes the average value of c™ in x and y directions, i.e. 


(*)?- l,fc + (*)& 


(<*)&-! + ( c ‘)^ 


\ c iL)hk — 0 i \ c iR)j,k — 


4. Numerical results 

For the sake of simplicity, we only give examples for the Poisson-Nernst-Planck system with two 
ions, i.e. ( 1 ). 

From the numerical scheme (8)-(10), it is easy to check that the truncation error of ( 8 ) and (9) are 
0(At + (Ax) 2 + (Ay) 2 ), while (10) has the truncation error 0((Ax) 2 + (Ay) 2 ). Thus the numerical 
scheme is expected to convergent with first order in time and second order in space. 

Example 1. Since it is not possible to find the exact solutions for the equations (1), we are now 
use the following argumented equations with exact solutions as a test problem: 

Pt = V ■ (Vp + pV0) + A, in n T = [0, T\ x Cl, (40) 

n t = V • (Vn - nV0) + f 2 , in Cl T = [0, T} x Cl, (41) 

—A <j> = p — n + c, in CIt = [0, T] x Cl, (42) 

where 

p = ( 3x 2 — 2a; 3 + 3 y 2 — 2 y 3 )e~ t , n = (a; 2 (l — a;) 2 + y 2 (l — y) 2 )e~ t , cj> = x 2 ( 1 — a;) 2 y 2 (l — y) 2 e _t 

are the exact solutions of (40)-(42), which satisfy the zero Neumann boundary conditions (3). And 
/i, / 2 , c are known functions which are given according to these exact solutions. 

We do the discretization of equations (40)-(42) as follows: 


m+i _ m 
Pj,k l J j,k 

= A hP ™+ 1/2 + (4(< fc 4C* +1/2 ) + WACfc +1/2 )} 


(43) 

At 

m+l _ m 
U j,k U j,k 

= A, l n™ +1/2 - (4(n™ t 4C +1/2 ) + E(</Ao;V ' N 

) + (h)7!k, 

(44) 

At 


-A,Cfc +1 = PTk 1 - k 1 + c#- 1 . 


(45) 


In the numerical computation, since the electric potential <f> is not unique up to a constant, we set 
electric potential at the first point to be the exact value at each time step in order to get unique 
solutions. It is easy to check that the truncation error of the above discretization scheme for the 
system (42) has the same order as the truncation error of the discretization scheme (8)-(10) for the 
problem (1). For this example, when numerically implementing of (43)-(45), we set (f> m+1 at one 
boundary point to be the exact value at each time step so that (ft m+1 is uniquely determined. 

Now we carry out the numerical convergence study for both space and time using (43)-(45). For 
spatial convergence, we set At = 0 . 000002 , and use 4 different spatial meshes Ax = Ay = 2Q 3 2 „ A h, 


9 























n = 0, • • ■ ,3, the final time is set to be T = 1.0. When At is sufficiently small, we compute the 
spatial convergence order according to 


order 1 = log 2 


\\uh(-,-,T) -u exact {-,-,T)\\ 

|| Uh(-, •, T) - U eX act (*5 ',T )\\’ 


(46) 


where Uh(-,-,T) is the numerical solution at time t = T using mesh h, u e xact {-, • ,T) is the exact 
solution at time t = T, and || • || is the spatial discrete norm. Table 1 shows the mesh refinement 
analysis for p , n, <j) using two different norms. One can see, the errors are decreasing when spatial mesh 
is refined, and it is second order convergent for both norms, which is expected from the truncation 
error analysis. 


Table 1 

Spatial mesli refinement analysis for the Poisson-Nernst-Planck system with zero Neumann boundary conditions 
(Cp = Ph Pexacti 6n — TLfi Tlexact ? — 4>h $ exact 5 At = 0.000002). 


h 11 11 2 or derl ll e plloo orderl ||e n || 2 orderl llenlloo orderl ||e</>|| 2 orderl || e ^>|| orderl 


1/20 

3.48e-4 

- 

7.65e-4 

- 

3.16e-3 

- 

1/40 

8.71e-5 

2.00 

1.94e-4 

1.98 

7.89e-4 

2.00 

1/80 

2.17e-5 

2.00 

4.91e-5 

1.98 

1.97e-4 

2.00 

1/160 

5.42e-6 

2.00 

1.25e-5 

1.97 

4.94e-5 

2.00 


3.24e-3 

- 

5.19e-3 

- 

5.74e-3 

- 

8.09e-4 

2.00 

1.63e-3 

1.67 

1.77e-3 

1.70 

2.02e-4 

2.00 

4.87e-4 

1.74 

5.23e-4 

1.76 

5.06e-5 

2.00 

1.40e-4 

1.80 

1.49e-4 

1.81 


Table 2 

Temporal mesh refinement analysis with h = 1/640 for the Poisson-Nernst-Planck system with zero Neumann bound¬ 
ary conditions. 


At ||ep || 2 order2 H^pIIqo order2 || || 2 or der2 llenlloo order2 || e <£|| 2 order2 \\ e <l>\\ order2 


1/40 

8.26e-3 

- 

1.31e-2 

- 

5.46e-4 

- 

8.61e-4 

- 

6.51e-3 

- 

7.40e-3 

- 

1/80 

4.14e-3 

1.00 

6.58e-3 

1.00 

2.75e-4 

0.99 

4.33e-4 

0.99 

4.05e-3 

0.69 

4.52e-3 

0.71 

1/160 

2.07e-3 

1.00 

3.29e-3 

1.00 

1.39e-4 

0.98 

2.18e-4 

0.99 

2.42e-3 

0.74 

2.67e-3 

0.76 

1/320 

1.04e-3 

1.00 

1.65e-3 

1.00 

7.11e-5 

0.97 

l.lle-4 

0.98 

1.41e-3 

0.78 

1.53e-3 

0.80 


For time convergence, we set Aa; = Ay = h = 1/640, and use 4 different time steps At = 40 ^ 2 „ , 
n = 0, ••• ,3, the final time is set to be T = 1.0. When h is sufficiently small, we compute the 
temporal convergence order according to 


order2 = log 2 


| -F) U eX act (*7 F) 11 

\\uAt/ 2 (-,-,T) - U exac t(-,-,T)\\’ 


(47) 


where UAt(-,-,T) is the numerical solution at time t = T using the time step At. Table 2 shows 
the time step refinement analysis for p,n,(j) using two different norms. One can see, the errors are 
decreasing when time step is refined, and it is first order convergent in time, which also is expected 
from the truncation error analysis. 

Example 2. We consider the equations (1) in the domain = [0,1] x [0,1] with zero Neumann 
boundary conditions (3), and initial conditions 


P(0 ,x,y) 
n(0,x,y) 



23 

144’ 


where the initial conditions are set to satisfy the conditions (2)-(4). 
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(a) electric potential energy E w.r.t. time 

0.1667 -,-T-T-T- 

0.1667 ■ 

0.1667 ■ 

g 0.1667 ■ 

2 

0.1667 ■ 

0.1667 ■ 

0 0.2 0.4 0.6 0.8 1 

(c) mass of n w.r.t. time 



(e) relative mass error of n w.r.t. time 





(b) mass of p w.r.t. time 





(d) relative mass error of p w.r.t. time 



(f) minimum p in Q w.r.t. time 



time 

(g) minimum n in £1 w.r.t. time 


Fig. 1. Numerical results for the Poisson-Nernst-Planck system with zero Neumann boundary conditions. 
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We carry out numerical computation with h = A t = 0.01 ,T = 1 using the scheme (8)- (10). 
In the numerical computation, we set electric potential at the first point to be zero at each time step 
in order to get unique solutions. Figure 1 gives the evolution of electric potential energy, mass of p , 
mass of n, the relative mass error of p, the relative mass error of n, minn(p), minn(n), respectively. 
One can see that the electric potential energy decays, and the mass of p, n is exactly conserved. 
Moreover, p and n always keep positive. All these results are consistent with the analysis in above 
section. 

5. Conclusions and discussions 

Prohl et al. [1] first proposed a fully implicit finite element method for the Poisson-Nernst-Planck 
system. Numerically, in order to get the rigorous mass conservation and electric potential energy 
decay properties, a fixed iteration method is needed for the fully implicit finite element scheme. In 
this paper, we develop a simple semi-implicit finite difference method for the Poisson-Nernst-Planck 
system, which can also preserve mass and electric potential energy identities numerically. The current 
method only needs to solve a linear system at each time step, which can be done very efficiently since 
the all coefficient matrices are symmetric banded. Furthermore, mesh refinement analysis shows that 
the method is second order convergent in space and first order convergent in time. And the method 
can be easily extended to the case of multi-ions. 

Since the Poisson-Nernst-Planck system is a nonlinear system, theoretical convergence analysis 
for the proposed numerical method will be a challenge task, we leave it as the future work. More¬ 
over, constructing simple and efficient numerical method which can preserve the entropy law of the 
Poisson-Nernst-Planck system is another future goal. 
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